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Abstract 

The  three-dirnensional  developrnent  of  flow  transition  in  botli  subsonic  and  supersonic 
Joukowsky  airfoil  boundary  layers  are  studied  by  direct  numerical  simulation  (DNS).  The 
r.umerical  simulation  is  performed  by  a  spatial  approach.  A  full  compressible  Navier-Stokes 
system  in  curvilinear  coordinates  is  employed  so  that  we  can  simulate  the  transition  around 
general  geometric  configurations.  The  numerical  results  agree  very  well  with  the  lineair  sta¬ 
bility  theory  (LST)  at  the  linear  growth  stage  for  both  primary  and  second  modes  in  the 
fiat  plate  boundary  layers.  The  whole  process  of  controlled  flow  transition  induced  by  blow- 
i.ig/suction  around  airfoils  is  simulated  by  directly  solving  the  N-S  system  with  Reynolds 
number  around  10®.  Some  differences  are  found  in  comparison  to  the  incompressible  coun¬ 
terpart,  and  some  new  phenomena  for  the  transition  around  airfoils  are  observed  which  at 
least  qualitatively  agree  with  physics. 
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1  Introduction 

In  recent  years,  there  has  been  a  continuous  and  rapid  improvement  in  understanding 
the  process  of  transition  from  laminar  to  turbulent  flow  in  incompressible  boundary  layers. 
Unlike  the  relatively  comprehensive  picture  of  transition  scenarios  in  incompressible  flows, 
nonlinear  effects  responsible  for  transition  at  high  speeds  are  still  very  much  a  mystery 
(Erlebacher  &  Hussaini,  1990).  In  the  absence  of  high-speed  experimental  data,  direct 
numerical  simulation  (DNS)  now  plays  an  important  role  in  understanding  high-speed  flow 
transition. 

Several  works  on  subsonic  boundary  layer  transition  have  been  published  recently  by 
other  researchers  (e.g.,  Fasel,  1990;  Thumm  et.  al.,  1990;  Masad  k  Nayfeh,  1992;  Masad 
k  Iyer,  1993;  Masad  k  Nayfeh,  1993).  As  the  power  of  today’s  computer  increases,  direct 
numerical  simulation  (DNS)  for  realistic  flow  at  certain  Reynolds  numbers  becomes  feasible. 
It  is  generally  agreed  that  compressibility  does  not  change  the  fundamental  physics  for 
strcarnwise  instabilities  at  subsonic  Mach  number,  but  it  still  has  many  differences  from 
incompressible  instability,  and  the  numerical  simulation  work  is  still  very  limited,  especially 
for  the  complex  geometries. 

For  supersonic  flow  transition,  DNSs  of  flat-plate  boundary  layers  over  a  wide  range  of 
Mach  numbers  (up  to  Mach  8  or  higher)  have  resulted  in  some  encouraging  quantitative 
agreements  (to  several  digits)  with  theories  (Pruett  et  al,  1995).  Several  researchers  have 
done  great  work  on  DNS  for  supersonic  boundary  layer  transition  (e.g.,  Erlebacher  k 
Hussaini,  1990;  Pruett  k  Zang,  1992;  Ei  Jler  k  Bestek,  1993;  Fasel  et  al,  1993;  Pruett 
et  al,  1995;  Adams  k  Kleiser.  1996;  Gno  et  al;  1996).  However,  the  high-speed  transition 
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is  very  sensitive  to  subtle  changes  in  inflow  disturbances  and  basic  state  flow  as  well.  The 
numerical  simulation  even  for  controlled  transition  in  general  geometries,  such  as  airfoils, 
still  remains  a  grand  challenge  up  to  date. 

An  alternative  approach  for  studying  flow  transition  is  the  PSE  (parabolized  stability 
equation)  method  (Herbert  &  Bertolotti,  1987;  Bertolotti,  1991).  This  approach  overcomes 
some  severe  limitations  of  the  traditional  linear  stability  theory  (LST)  by  accounting  not 
only  for  streamwise  variations  but  also  for  inhomogeneous  and  nonlinear  terms.  Unfortu¬ 
nately,  the  cost  of  PSE  increases  much  faster  than  that  of  DNS  as  the  number  of  spanwise 
modes  is  increased.  PSE  thus  is  limited  in  practice  to  investigations  of  narrow-band  forcing 
(Pruett  et  al,  1995).  Furthermore,  the  parabolization  limits  information  about  downstream 
feedback,  which  is  physically  inaccurate.  Other  methods,  like  compressible  linear  stability 
theory  (LST)  (Mach,  1984)  and  secondary  instability  theory  (SIT)  (Ng  &  Erlebacher,  1992) 
can  also  provide  limited  information  for  flow  transition,  but  are  far  away  from  practical 
application. 

Therefore,  in  order  to  really  solve  the  practical  problems,  an  efficient  and  accurate  solver 
for  the  full  compressible  time-dependent  Navier-Stokes  system  must  be  developed.  To  date, 
Navier-Stokes  codes  fall  roughly  into  two  classes  depending  upon  the  application  (Pruett  et 
al..  1995),  aerodynamic  codes  which  can  handle  complex  geometry  but  are  of  low  accuracy, 
and  DNS  codes  which  are  of  high  accuracy  but  can  only  handle  simple  geometries,  like 
channel,  flat  plate,  cind  cone.  Our  goal  is  to  develop  a  DNS  code  which  is  of  high  accuracy 
and  can  also  handle  complex  geometry,  so  that  we  can  simulate  the  real  time-dependent 
compressible  transitional  and  turbulent  flow  around  aircraft  components  or  subsystems, 
and  eventually  Air  Force  flight  vehicles  at  all  speeds. 

Liu  et  al.  (1996a,b)  developed  an  approach  which  enables  us  to  do  the  simulation 
of  the  whole  process  of  transition  in  the  incompressible  boundary  layers  for  airfoils.  For 
the  compressible  flow,  we  can  still  adopt  some  of  the  techniques  from  our  previous  work. 
Though  the  implicit  schemes  have  better  numerical  stability,  they  consume  much  more 
computer  resources  than  explicit  schemes,  especially  for  complex  geometry.  For  unsteady 
problems  the  time-step  is  very  small,  so  that  explicit  schemes  can  also  work  well.  In 
addition,  explicit  code  is  easy  to  vectorize  and  parallelize.  In  this  study,  we  use  explicit 
methods. 

For  subsonic  transition  cases,  an  explicit  time-steeping  code  is  generated  using  a  sixth- 
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order  central  difference  for  the  convection  terms,  a  fourth-order  central  difference  for  the 
viscous  terms,  and  a  four  stage  fourth-order  Runge-Kutta  method  for  time  integration. 
Under  this  frsnne,  the  CFL  number  can  be  increased  to  around  2.8,  so  that  the  time- 
step  size  for  numerical  stability  is  not  too  restricted.  For  supersonic  transition  cases,  we 
generate  our  code  for  explicit  time-stepping  by  adopting  a  sixth-order  compact  central 
difference  (Lele,  1992)  for  both  convection  and  diffusion  terms,  and  a  three-stage  third- 
order  compact  Runge-Kutta  method  (Wray,  1986)  for  time  integration. 

The  accuracy  of  the  code  is  examined  at  the  linear  instability  stage  by  comparing  the 
numerical  results  with  the  linear  stability  theory  (LST)  at  Moo  =  0.5  and  Moo  =  4.5.  The 
perturbation  is  introduced  by  wall  suction/blowing  upstream.  After  a  certain  distance  of 
development,  the  least-stable  mode  is  picked  up  and  behaves  the  same  ss  shown  by  LST. 
Two  typical  cases  of  the  whole  process  of  transition  in  the  subsonic  flat  plate  boundary 
layers,  K-type  cind  H-type  breakdown,  are  calculated  to  ensure  the  robustness  of  the  code. 
Transition  over  2-D  Joukowsky  airfoils  is  then  numerically  simulated  for  both  subsonic  and 
low  supersonic  Mach  numbers.  Some  new  phenomena  are  observed. 


2  Governing  System 


The  3-D,  compressible,  time-dependent  Navier-Stokes  equations  are  considered  as  the 
governing  system  in  this  work.  On  the  general  curvilinear  coordinate  system  (^,  77,  C),  these 
equations  can  be  expressed  as  follows: 
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Here.  V=  (n,v,w)'^.  The  four  dimensionless  parameters  resulting  from  the  nondimension- 
alization  are  defined  as: 


4 


Re  = 
Pr  = 


Poo^oo^ 

f^oo 

CpPoo 


7  =  TT 


P 

a’ 


(3) 


where  Rg,  Cp,  and  Cv  are  the  ideal  gas  constant  and  the  specific  heats  at  constant  pressure 
and  constant  volume,  respectively.  Generally,  we  choose  Pr  =  0.72  cind  7  =  1.4  for 
subsonic  and  low  Mach  number  supersonic  flows.  Viscosity  is  normalized  by  the  viscositv 
at  freestream,  and  is  assumed  to  vary  according  to  Sutherland’s  law.  This  yields 


1+C 
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(4) 


The  lengths  are  scaled  by  some  reference  length  /*.  For  the  flat  plate  boundary  layer  flow, 
we  use  the  boundary  layer  displacement  thickness  at  the  specified  location,  where  the 
disturbance  is  introduced.  In  general,  with  the  dimensionless  form,  8*  is  defined  as 


S”  =  f  {l-pu)dy. 

JQ 

For  the  airfoil  boundary  layer  flow,  the  half-thickness  of  the  airfoil,  h, 
length. 


(5) 

is  used  as  the  reference 


lo  reduce  the  machine’s  round-off  error,  we  decompose  the  total  flow  into  two  parts: 
steady  base  flow  and  unsteady  perturbation  flow  for  the  current  work  (in  general,  the  base 
flow  can  also  be  unsteady  flow).  Assume 
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where  the  subscript  “0”  denotes  the  basic  steady  state  flow,  and  prime  (')  denotes  the 
perturbation.  The  resulting  governing  system  in  the  perturbation  form  is: 
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In  the  above  governing  system,  the  vectors,  the  shear-stress  tensor,  and  the  heat-flux  vector 
are  given  as  follows: 
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Again,  the  subscript  “0”  denotes  the  basic  steady  state  flow,  prime  denotes  the  perturba¬ 
tion,  V'  =  (u',u',u;')^,  and  Vq  =  («o.Vo,ti^o)^- 

With  solved,  the  perturbation  components  can  be  retrieved  from 
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3  Numerical  Method 


Because  of  the  sensitivity  of  physics,  DNS  codes  must  meet  the  requirements  of  both 
high  accuracy  and  strong  numerical  stability.  Many  high-order  discretizations  have  been 
developed  by  other  researchers,  for  example,  spectral  methods  (Orszag,  1971),  compact 
finite  difference  methods  (Lele,  1992),  upwind-biased  differences  (Rai  &:  Moin,  1993),  and 
standeird  high-order  central  differences.  Because  the  boundary  treatment  for  the  compact 
JifTeieace  scheme  is  sensitive,  we  use,  for  the  subsonic  cases,  the  standard  sixth-order 
central  difference  for  the  convective  terms,  fourth-order  central  difference  for  the  viscous 
terms.  In  the  supersonic  cases,  we  choose  the  sixth-order  compact  finite  difference  method 
(Lele,  1992)  for  both  convection  and  diffusion  terms. 

For  time-integration,  most  people  use  the  so-called  compact  low-storage  third-order 
Runge-Kutta  scheme  (Wray,  1986)  for  simulating  the  transition  problem.  This  method  is 
efficient  for  supersonic  flow  because  with  the  same  CFL  number,  the  real  time-step  size  in 
supersonic  cases  is  relatively  larger  than  in  the  subsonic  situation.  So  in  the  supersonic 
cases,  we  adopt  this  method.  For  the  subsonic  cases,  we  adopt  the  classic  four-stage 
fourth-order  Runge-Kutta  scheme  with  sixth-order  central  difference  in  space  directions. 
By  carefully  rearranging  the  arrays  in  the  code,  we  can  arrange  that  this  fourth-order 
Runge-Kutta  scheme  also  has  low  storage. 


We  can  rewrite  the  system  into  the  semi-discrete  form: 


-U:,.,-fR(U',.,)  =  0, 


where  R(U'  _,j;.)  is  the  function  defined  by 


(10) 


R(u:j.,)  =  [Lc  +  LzjjU',,.  „ 


(11) 


where  Lc  and  Ld  are  discretized  operators  for  convection  and  diffusion,  respectively.  The 
fourth-order  Runge-Kutta  scheme  used  to  solve  a  system  of  ODE  corresponding  to  the 
discretized  Navier-Stokes  equations  can  be  expressed  as  follows: 


(12) 


U'(3)  _  u'(°)  - 

U'(4)  = 

6 

U'(n+1)  _ 


Here,  =  R(U'^^^),  the  superscript  n  denotes  the  time  level  at  nAt,  and  the  mesh 

indices  {i,j.k)  associated  with  the  solution  vector  and  time  step  At  are  suppressed  for 
convenience. 


4  Boundary  Conditions 

The  above  system  can  be  solved  with  proper  specification  of  boundary  conditions.  In 
general,  all  flow  quantities  can  be  specified  at  the  inflow  boundary.  The  prescription  of  the 
inflow  boundary  conditions  depends  on  the  way  that  the  disturbances  are  introduced  into 
the  bcise  flow.  We  can  use  either  LST  or  PSE  to  provide  the  eigenvectors  at  the  inflow. 
However,  we  think  a  periodic  blowing/suction  at  the  wall  is  more  natural. 

Figure  1  gives  a  schematic  description  of  the  computational  domain.  Supposing  the 
blowing/suction  strip  is  located  relatively  far  away  from  the  inflow  boundary,  we  can  assume 


f'{xo,y,z,t)=^0, 


(1.3) 


where  /'  stands  for  the  variables  u',v',w',p',  and  p\  respectively. 

For  the  boundary  conditions  at  the  wall,  we  use  no-slip  and  isothermal  conditions; 


“'(x,t/waH. -.0  =  w'{x,y^alhZ,t)  =  T'{x ,  y^alh  Z ,  t)  =  0. 


(14) 


The  wall  normal  velocity  component  is  prescribed  cis  a  function  of  x,z,  and  t  to  generate 
the  disturbances  in  the  computational  domain  (Rist  &  Fasel,  1995); 


v{x,y^aii,z,t) 


0  X  <  ii  or  X  >  X2. 


In  this  work,  we  use 


Kiz:,z,t) 


+ 


£2(f  sin 
€3d  sin 


( 27r(x  -  xi)\ 

V  (a:2-xi)  / 

/27r(x  -Xi)\ 

V  (x2-x,)  J 


sin{u>2di) 

cos(P  z)sin{u)3dt), 
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where  e^d  and  e^d  are  the  amplitudes  for  2D  and  3D  disturbances,  respectively;  (5  is  the 
controlled  spanwise  wave  number;  and  uj2d  and  uzd  are  the  frequencies  for  2D  and  3D 
disturbances.  This  disturbance  represents  a  combination  of  one  2D  wave  and  a  pair  of 
oblique  waves. 

In  the  spanwise  direction,  a  periodic  boundary  condition  is  imposed. 

In  the  far  field,  an  exponential  decay  condition  (Fasel  et  al,  1990)  together  with  the 
sponge  layer  treatment  (Collis  Lele,  1996)  are  employed  to  satisfy  the  non-reflection 
requirement.  The  exponential  decay  condition  yields 

f  =  f  _  (15) 

Here, 

ys;  ’ 

where  q,,  is  the  real  part  of  the  streamwise  wave  number. 

The  sponge  layers  are  also  embedded  in  the  inflow  and  outflow  sections  (see  Figure  1) 
to  eliminate  the  reflection  in  the  streamwise  direction.  For  subsonic  flow,  the  inflow  sponge 
is  necessary,  while  for  supersonic  flow,  the  inflow  sponge  can  be  removed. 

For  the  Runge-Kutta  type  method,  the  sponge  is  loaded  at  each  single  stage,  thus,  (12) 
needs  to  be  modified  by  adding  an  exponentially  decaying  source  term  in  space  directions 
to  the  right  hand  side: 

^U'(x,  y,  z)  -f  R(U'(x,  y,  z))  =  -/d(x,  y)U'(x,  y,  z).  (16) 

This  artificial  term  generates  a  so-called  stiff  differential  equation  which  needs  to  be  treated 
implicitly.  This  yields: 

U'{o)  ^  u'(n)^ 

U'C)  =  lu'W-I^Rl'll/ll+l^/^x,!,)], 

U'W  =  |u'm-^R»il/li  +  y/4i,!/)l. 

U'(3)  ^  +  Atfd{x,y)], 

U'(4)  =  [U'(0)  _  —{-rW  +  -f  2R<^^  R<^>)]/[1  +  At /^(x,  y)], 

6 

U'("+i)  =  (17) 
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The  sponge  function  is  given 

/  C  /  \  A  f 


fd{x,  y)  = 


/d3(x,y)  = 


Mx,y)  = 
/^s(a:,y)  = 
0 


X  e  [xo,x^i), 

X  €  (Xj2iXmaa;]) 

X  €  [x,i,x,2], 

X  €  [xqj  Xji), 

X  6  (x«2)Xma«]) 

elsewhere 


y  e  [o,y.) 
y  €  [0,y,) 
y  €  (y»  j  ymox] 
y  e  (y„ymar] 
y  ^  (y*»  ymax] 


Here,  A^i,  As2,  A^s,  Ngj,  N^z,  and  Nss  are  prescribed  constants.  In  this  work,  we  choose 
.4si  =  As2  =  Aj3  =  10,  and  iV,i  =  N^z  =  Nss  —  3. 


5  Grid  mapping 

To  assure  that  the  numerical  schemes  of  DNS  have  high  accuracy,  we  need  to  make  the 
accuracy  of  grid  transformation  schemes  two  orders  higher  than  that  of  the  flow  scheme 
itself.  In  this  paper,  we  only  use  a  flat  plate  or  a  Joukowsky  airfoil  with  zero  attack  angle  as 
our  investigation  targets,  so  we  are  able  to  use  analytical  mapping  to  transfer  the  physical 
domain  to  the  computational  domain. 

In  subsonic  and  low  supersonic  Mach  number  flows,  the  primary  mode  dominates  the 
transition  process,  while  as  Mach  number  goes  to  over  Moo  =  3,  the  second  mode  (Mach 
mode)  becomes  dominant.  This  requires  use  of  different  types  of  grid  structure  for  the 
above  two  kinds  of  flows. 


5.1  Grids  for  primary  instability 


For  the  problem  of  the  flat  plate  boundary  layer  transition,  the  grid  mapping  function 
we  used  is 

r]mazy{f^  +  ymax)  I  ^  in 

-\l - • 

X 


vi^^y)  = 


(19) 


ymaxicr  +  yyjxin/x) 

Here,  x,n  is  the  location  where  Re*  is  defined.  Usually,  we  choose  =  X2  (see  Figure 
1).  For  the  parallel  bcise  flow,  the  boundary  layer  thickness  does  not  grow,  so  the  above 
transformation  becomes 

^maxy{^  “h  ymax) 


y{y)  = 


(20) 


Vmaxicr  +  y) 

.■\  typical  grid  for  2D  flat  plate  boundary  layer  transition  is  shown  in  Figure  2  (a).  Note 
that  we  use  straight  lines  (20)  for  x  <  .r„,.  and  parabolas  (19)  for  x  >  x,„.  This  avoided 
using  too  small  time-steps  required  for  the  CFL  condition. 


(18) 
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5.2  Grids  for  second-mode  instability 

For  high-speed  (Moo  ^  3)  flow  transition,  high  gradients  are  usually  concentrated  near 
the  wall  and  near  the  critical  layer.  With  the  aid  of  LST,  we  can  estimate  the  location  of 
critical  layers,  and  then  use  an  analytical  mapping  which  clusters  grid  points  in  both  wall 
and  critical  layer  regions.  For  the  flat  plate  with  parallel  base  flow,  we  use  the  following 
mapping,  which  is  similar  to  that  used  by  Pruett  et  al.  (1995). 

1.  In  computational  space: 


2.  Map  T]  to  [-1,1]: 


where  r/1/2  =  rj^j/2. 


m  =  (i  - 1), 

i  =  l,2,---,iVJ. 


-  rg  -  ?7i/2 
^1/2 


(21) 


(22) 


3.  There  are  five  parameters  for  the  mapping:  t/OTai,!/i/2,J/o,  Ay,  and  r.  As  shown  in 
Figure  2(b),  the  parameter  ymax  is  the  height  of  the  computational  domain;  yi/2 
defines  a  point  such  that  half  of  the  total  grid  points  lie  between  the  wall  eind  this 
point;  yo  is  usually  chosen  near  the  critical  layer,  and  Ay  and  r  are  used  to  adjust 
the  stretch  strength  near  yo.  The  next  step  maps  [-1,1]  onto  itself  but  performs  a 
clustering  of  grids  about  yo; 


4.  Exponential  transformation  maps  [-1,1]  onto  a  physical  interval  [0,y„,a 


(23) 


yW  =  yma 


-1  <  V'  <  1, 


a2-  1 

ymax  "t"  yjy'max  4yi/2(ymax  yi/2) 
2yi/2 

The  remaining  unknown  quantities  in  (23)  are  determined  ais  follows: 


a  = 


(24) 


;f^(a^-l))/log(a)-l, 

ymax 

Arp  =  Ayol{dy/drp), 

Vm  zx 


rpo  =  log(l  + 
A0  = 
dy/drp  = 


A4'o+l) 


(a^-l) 


log(a). 
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Equation  (23)  is  solved  by  Newton’s  method  to  obtain  ‘ip,  and  then  (24)  is  used  to  obtain  the 
grid  coordinates  in  the  physical  interval  [Ojt/maij-  For  the  nonparallel  flat  plate  boundary 
layer,  the  local  boundary  thickness  will  be  changed  according  to  ^  =  Sinyjxjxin,  so  we  let 

y  =  Vinsfxj^ 


where  t/,n  is  the  normal-direction  coordinate  we  obtained  above  (at  inflow).  All  the  Jacobian 
metrics  are  obtained  analytically. 

5.3  Conformal  Mapping  for  the  Joukowsky  Airfoil 

For  the  airfoil  boundary  layer  flow,  we  now  restrict  our  problem  to  be  zero-degree 
attack  angle.  For  Joukowsky  airfoils,  conformal  mapping  is  used  to  obtain  the  grids  and 
transformation  metrics  analytically  (Liu  et  al.,  1996  a,b). 

The  grid  generation  process  can  be  described  as  follows  (Figure  3): 

1-  {(yrj)  {(Q,r]o): 

fo  = 

(^A 

T)o  =  TJ  . 

This  step  performs  a  uniform  stretch  in  the  streamwise  direction  so  that  the  length 


of  solid  wall  is  2a.  This  implies 


(^A  = 


6. 

2a' 


2-  (^o^t/o)  (^i)^i)- 


-a, 


6  = 

(^CVmax  +  Vl  max  iv  max  -voY 


(Tb  +a  -  l^ol 

<^CT?lmox^0 


The  first  formula  performs  a  clustering  of  grids  about  =  —a,  and  the  second  one 
performs  a  clustering  of  grids  about  solid  wall  {rji  =  0). 


3-  (6i^i)  (a:i,yi): 

By  using  the  Joukowsky  transformation 

6  =  ^(1  + 


m 


yi 

2  ^  xl  +  yl 
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we  can  evaluate  xi  and  yi. 


4.  (a:2,y2): 


Xi 

Vi 


Xi  - 


k'a 

F+T’ 


yi- 


This  step  performs  a  shift  of  j/i  axis  in  the  streamwise  direction. 


5.  (x,y): 


X 

y 


CTEXiil  + 


<^£y2(i  - 


{k'  +  1)2  xl  +  yl^' 

— f! _ 1_) 


This  step  will  generate  the  grids  we  required.  Here,  A:'  <  1  is  used  to  change  the  ratio 
of  half-thickness/chord-length.  For  example,  with  k'  =  0.0116,  the  ratio=l/16.  .Also, 
we  use  cte  to  normalize  the  airfoil  and  make  the  half-thickness  equal  to  1. 


With  the  grid  generated,  we  can  also  obtain  the  Jacobian  metrics  between  the  (^, /?) 
plane  and  the  (x,  y)  plane  analytically. 


6  Filtering 

Because  of  the  use  of  the  central  difference  scheme,  the  two-point  “saw-tooth”  oscilla¬ 
tions  will  generally  be  generated,  especially  in  the  low  viscous  region.  Those  oscillations  can 
induce  some  spurious  physical  waves.  To  avoid  this  phenomenon,  we  can  either  put  some 
scheme  viscosity  through  the  biased  upwinding  scheme  (Rai  &  Moin,  1993),  or  through 
the  two-step  MacMormack-like  dissipative  method  (Gottlieb  k  Turkel,  1976),  or  through 
explicit  filtering.  Due  to  the  simplicity  of  explicit  filtering,  we  adopt  it  in  our  code.  A 
sixth-order  filter  is  applied  every  10  time-steps  in  the  wall-normal,  direction,  and  every  20 
steps  in  the  streamwise  direction  for  the  subsonic  caises.  During  the  study  of  supersonic 
transition,  we  find  that  only  the  wall-normal  direction  filtering  is  necessary.  A  compact 
sixth-order  filtering  is  used  every  10  time-steps  in  that  direction. 
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7  Intergrid  dissipation 


Current  supercomputers  are  still  far  away  from  resolving  all  the  length  scales  for  fully 
developed  turbulent  flow,  while  the  use  of  models  to  represent  those  unresolved  scales 
can  always  cause  the  case-dependent  problem.  But,  without  using  any  model,  most  non- 
dissipative  DNS  schemes  will  blow  up  because  of  the  Itick  of  enough  dissipative  mechanism 
of  the  codes.  The  “intergrid  dissipation”  developed  by  Liu  &  Liu  (1995)  plays  a  positive 
role  in  our  code.  By  using  the  difference  of  flow  quantities  between  two  grid  levels,  we  are 
now  able  to  eliminate  those  scales  badly  described  by  the  current  grid,  and  provide  some 
dissipation  caused  by  those  subgrid  scales. 

A  brief  description  is  as  follows.  Suppose  is  the  quantity  obtained  on  the  global  grid 
G^.  We  first  restrict  it  to  the  next  coarse  level  grid  and  then  obtain  the  new  value 
of  on  the  global  grid  G^.  However,  if  we  simply  adopt  this  the  effects  of 

small  scales  will  all  be  ignored.  Therefore,  a  dynamic  weight  function,  cr,  is  introduced  to 
maintain  the  effect  of  small  eddies.  This  yields 

fh{new)  ^  ^jh^j2hjh{old) 

a  =  min{2/c,  1}, 

«  =  ^(Po  +  1p'1)(«'^  + 

Here.  is  a  full- weighting  restriction  operator,  cind  is  a  bi-linear  interpolation  operator. 
For  the  detail  of  intergrid  dissipation,  see  Liu  Liu  (1995). 

8  Computational  Results  and  Discussion 

8.1  Subsonic  transition 

8.1.1  Flat  plate 

The  code  is  first  validated  for  the  Mach  0.5  flat  plate  boundary  layer  flows.  The  base 
flow  is  obtained  by  solving  the  compressible  similarity  system  (Stewartson,  1964)  using  a 
shooting  method.  The  solid  is  assumed  to  be  adiabatic  for  the  base  flow  and  isothermal  for 
the  perturbation.  The  Reynolds  number  is  set  to  Re’  =  875  (based  on  <J"),  and  u>2d  =  0.1. 
Compressible  linear  stability  theory  provides  an  eigenvalue  a  =  0.2636  —  iO. 005623. 

Tlie  whole  computational  domain  is  set  to  15  T-S  wavelengths  to  ensure  the  development 
of  the  least  stable  mode.  The  grid  is  16/wave  x51  x  1,  with  the  last  wavelength  used  as 
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the  buffer  domain  (sponge  layer),  and  the  stretch  parameter  in  (19)  is  cr  =  2.5.  Also  one 
T-S  period  is  divided  into  1000  time  steps.  The  base  flow  is  assumed  to  be  parallel.  Figure 
4(a)  depicts  the  disturb2mce  amplitude  of  u'  and  v',  showing  that  the  least  stable  mode  is 
picked  up  very  well  and  grows  as  LST  predicted.  Also,  we  can  see  that  the  sponge  layers 
successfully  eliminate  the  reflecting  waves  in  all  directions  to  keep  the  physical  domain  very 
clean.  Eigenfunctions  of  this  case  are  also  compared  with  LST  and  found  to  agree  very 
well  (Figure  4(b)). 

Both  K-type  and  H-type  transitions  in  2-D  subsonic  flat  plate  boundary  layers  are  then 
simulated  by  using  the  same  code. 

For  the  K-type  transition,  we  use  a  193  x  41  x  32  grid,  which  includes  an  inflow  sponge 
(first  16  streamwise  grid  points),  a  suction/blowing  slot  (next  16  streamwise  grid  points), 
an  outflow  buffer  domain  (the  last  16  streamwise  grid  points),  and  an  approximately  9 
T-S  wavelengths  physical  domain.  The  stretch  parameter  in  (19)  is  set  to  <r  =  2.5.  The 
Mach  number  is  Moo  =  0.5,  and  Too  =  2SSK.  The  height  of  the  computational  domain 
is  30  (based  on  the  S*  at  the  end  point  of  suction/blowing  slot,  Xi).  The  amplitude  of 
disturbance  is  set  to  ejd  =  0.02,  esd  =  0.008,  the  controlled  spanwise  wave  number  is 
P  =  0.2.  and  the  angular  frequency  u  =  0.1.  One  forcing  period  is  divided  into  1000  time 
steps.  It  takes  around  8  Cray  C-90  CPU  hours  (one  processor)  to  run  the  code  for  10,000 
time  steps.  As  in  our  previous  work,  our  DNS  code  can  keep  running  without  blowing 
up,  and  the  process  of  breakdown  is  simulated  up  to  the  limit  of  current  grid  resolution. 
Figure  5  depicts  the  instantaneous  contour  plots  of  perturbation  streamwise  vorticity  (tu^) 
and  spanwise  vorticity  (a;,)  on  the  y"  =  0.18595,  0.48165,  0.85052  {x,z)  planes  after  10 
forcing  periods.  It  shows  that  the  transition  process  is  very  similar  to  that  of  incompressible 
flow,  except  that  the  “legs”  of  the  so-called  lambda  waves  are  longer.  Also,  even  after  the 
lambda  waves  break  into  smaller  scale  eddies,  the  major  splitting  appears  in  the  spanwise 
direction,  while  the  length  scale  in  the  streamwise  direction  is  relatively  longer. 

For  the  H-type  transition,  a  265  x  45  x  42  grid  is  employed.  In  the  streamwise  direction, 
the  whole  domain  is  divided  into  a  16  point  inflow  sponge,  a  24  point  suction/blowing  slot, 
a  24  point  outflow  buffer  domain,  and  a  physical  domain.  The  stretch  parameter  is  cr  =  2.2. 
The  Mach  number  is  set  to  Moo  =  0.8,  and  Too  =  28SK.  The  height  of  the  computational 
domain  is  25  (based  on  the  5’  at  the  end  point  of  suction/blowing  slot).  The  amplitude 
of  disturbance  is  set  to  S2d  =  0.012,  £3^  =  0.004,  the  controlled  spanwise  wave  number  is 
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(i  =  0.22,  and  the  angular  frequency  uj^d  =  0.15075.  One  forcing  period  is  divided  into 
550  time  steps.  It  takes  around  15  Cray  C-90  CPU  hours  (one  processor)  to  run  the  code 
for  12  forcing  periods.  Figure  6  shows  the  instantaneous  contour  plots  of  perturbation 
strccunwise  vorticity  (a;*)  and  spanwise  vorticity  (w,)  on  the  y*  =  0.1471,  0.3768,  0.6563 
(x,z)  planes  after  11  forcing  periods.  Basically,  the  transition  process  is  very  similar  to 
that  of  incompressible  flow,  with  staggered  structure  of  the  lambda  waves  at  the  secondary 
instability  stage.  Again,  the  “legs”  of  the  lambda  waves  are  relatively  longer.  One  difference 
from  the  incompressible  H-type  transition  is  that  the  staggered  lambda  waves  can  only 
be  observed  in  a  very  short  regime  for  the  compressible  case.  After  that,  the  staggered 
structure  is  contaminated  by  the  invading  lambda  waves  from  upstream,  and  a  branching 
structure  can  be  observed  clearly  from  the  vorticity  contours.  However,  the  structured 
pattern  can  be  observed  in  a  relatively  long  regime,  so  the  compressibility  has  an  overall 
effect  of  stabilizing  the  subhaxmonic  modes,  and  makes  the  transition  zone  longer. 

Since  the  perturbation  is  introduced  through  blowing/suction,  many  disturbance  modes 
may  coexist;  this  can  cause  our  simulations  to  be  more  or  less  like  the  detuned  type  (Ng  h 
Erlebacher,  1992). 

8.1.2  Joukowsky  airfoil 

The  first  case  here  is  the  fundamental  breakdown  of  a  subsonic  Joukowsky  airfoil 
boundary  layer  transition.  Note  that  now  the  reference  length  is  the  half-thickness  of  the 
airfoil,  h.  The  parameters  are  set  as  follows: 


Moo 

0.5, 

T„ 

= 

288A', 

Rth. 

= 

4000  (Ae^  =  160,000), 

Pr 

= 

0.72, 

LJ 

= 

0.8, 

P 

= 

1.3, 

^2d 

= 

0.015, 

£3d 

= 

0.01, 

a 

= 

0.25  , 
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The  chordlength/half-thickness  is  40,  which  makes  the  Reynolds  number  (based  on  the 
chordlength)  Rei  =  160,000.  A  281  x  41  x  31  grid  is  employed.  The  wing  is  placed  from 
streamwise  grid  point  41  to  261,  and  the  blowing/suction  slot  is  embedded  between  point 
71  and  point  79.  The  height  of  the  computational  domain  is  around  I6h,  and  the  crest 
point  of  the  airfoil  is  located  at  the  position  of  25%  chordlength.  Figure  7  gives  a  part  of 
the  grid  for  this  case.  For  the  time  stepping,  one  forcing  period  is  divided  into  1200  time 
steps. 

For  the  airfoil  boundary  layer  transition,  one  major  issue  is  that  it  contains  both  fa¬ 
vorable  and  adverse  pressure  gradients  in  the  streamwise  direction.  Before  the  disturbance 
reaches  the  crest  point  of  the  airfoil,  the  favorable  pressure  gradient  will  stabilize  the  flow, 
and  thus  the  amplitude  of  disturbance  is  damped.  After  the  disturbance  has  passed  the 
crest  point,  the  gradually  increased  adverse  pressure  gradient  will  destabilize  the  flow,  and 
the  amplitude  of  disturbance  is  expected  to  grow  faster  than  that  in  the  zero-pressure- 
gradient  flat  plate.  This  phenomenon  can  be  observed  clearly  in  the  2-D  case.  Figure 
8  depicts  the  maximum  perturbation  rms  amplitude  in  the  computational  domain.  The 
parameters  we  used  are  the  same  as  above,  except  e2d  =  0.001,  e3d  =  0.  Clearly,  the  pertur¬ 
bation  damps  much  faster  than  that  in  the  flat  plate  catse,  until  it  reaches  the  crest  point 
of  the  airfoil.  After  that,  a  very  strong  amplification  shows  up.  This  can  cause  breakdown 
when  3-D  disturbances  are  introduced. 

Figure  9  gives  the  contour  plots  of  perturbation  u\  v\w\  for  the  above  fun¬ 

damental  transition  case.  Similar  to  the  flat  plate  boundary  layer  case,  it  first  generates 
lambda  waves;  then  the  lambda  waves  break  into  smaller  scale  vortices,  and  cause  the  tran¬ 
sition.  One  difference  is  that,  in  the  present  case,  the  perturbation  u'  has  a  relatively  strong 
anisotropic  structure  (Figure  9  (a)).  Till  now,  we  have  not  found  an  answer  to  explain  this 
phenomenon.  Therefore,  a  more  detailed  investigation  needs  to  be  done  later. 

The  final  c^se  studied  for  subsonic  boundary  layer  transition  for  a  2-D  airfoil  is  the 
subharmonic  transition  case.  We  have  the  following  setups: 

Moo  =  0.6, 

ReH  =  13400  (/?eL  =  670,000). 


Pr 

= 

0.72, 

— 

1.206, 

^Zd 

= 

0.603, 

= 

1.76, 

^2d 

= 

0.01, 

^3d 

= 

0.0075, 

a 

= 

0.15, 

z 

€ 

A  361  X  45  X  41  grid  is  used,  in  which  the  wing  is  located  between  the  41st  and  321st 
streamwise  grid  point.  The  height  of  the  computational  domain  is  about  20  h,  and  the 
chordlength/half-thickness  is  50.  The  blowing/suction  slot  is  located  between  the  81st  and 
93rd  grid  points.  For  the  time  stepping,  one  forcing  period  (based  on  u;2d)  is  divided  into 
1000  time  steps.  It  takes  around  25  Cray  C-90  CPU  hours  (one  processor)  for  our  code 
to  run  for  11  forcing  periods.  The  contour  plots  of  perturbation  u',  v',  w',  p',  p'  for  this 
case  after  11  forcing  periods  are  given  in  Figure  10.  As  mentioned  by  Ng  &  Erlebacher 
(1992),  the  subharmonic  modes  are  the  most  unstable  modes.  This  results  in  earlier  onset 
of  transition.  In  fact,  the  staggered  structure  of  lambda  waves  is  much  weaker  than  that 
in  the  incompressible  flow.  Even  after  a  relatively  long  distance,  we  cam  still  find  that  u' 
and  p'  have  large  length  scale  in  the  streamwise  direction,  while  the  length  scale  in  the 
spanwise  direction  becomes  very  small.  This  could  cause  a  strong  anisotropic  structure  for 
the  turbulent  flow  around  airfoils. 

8.2  Supersonic  transition 

8.2.1  Code  validation 

The  code  for  supersonic  transition  first  is  validated  for  the  Mach  4.5  flat  plate  boundary 
layer  flows.  The  base  flow  is  obtained  by  solving  the  compressible  similarity  system  (Stew- 
artson,  1964)  using  a  shooting  method.  The  Reynolds  number  is  set  to  Re*  =  10,000, 
Pr  =  0.7,  Too  =  61.111,  and  wjj  =  2.0467.  The  constant  C  in  Sutherland’s  law  is 
C  =  110.33.  Compressible  linear  stability  theory  provides  an  eigenvalue  a  =  2.248— iO. 0308. 

A  16/wave  xl41  x  1  grid  is  used,  with  the  last  wavelength  used  as  the  buffer  domain. 
.Again,  one  T-S  period  is  divided  into  200  time  steps.  For  blowing/suction  case,  Ei/?ler  and 


19 


H.  Bestek  have  shown  that  a  second  mode  wave  and  an  additional  wave,  which  are  called 
“multiple- viscous-solution” ,  will  be  excited.  Since  the  “viscous-solution”  is  weakly  damped, 
the  second  mode  becomes  dominant  due  to  its  large  amplification  rate  if  the  computation 
domain  is  long  enough.  We  reach  this  goal  by  using  the  following  approach:  pick  up  all  the 
eigenfunctions  just  before  the  buffer  domain,  using  them  as  the  inflow  conditions,  instead 
of  using  the  blowing/shction  to  excite  the  perturbation  for  the  next  computation.  Figure 
11  depicts  the  disturbance  amplitude  of  u'  and  v',  showing  that  the  least  stable  mode  is 
picked  up  very  well  and  grows  almost  the  same  as  LST  predicted.  Though  we  still  can 
observe  a  long  wave  in  the  amplitude  curve,  it  is  very  weak. 

A  wall-normal  distribution  of  the  eigenfunctions  for  the  Mach  4.5  second  mode  is  also 
obtained  and  compared  at  Re*  =  10, 000  in  Figure  12.  It  can  be  observed  that  even  with 
blowing/suction,  our  code  can  still  pick  up  the  least  stable  mode  which  matches  the  result 
obtained  from  LST  very  well. 

8.2.2  Joukowsky  airfoil 

Though  the  second  mode  is  dominant  for  high  Mach  numbers,  the  transition  zone  is 
relatively  long.  On  the  other  hand,  for  low  supersonic  Mach  number,  the  viscous  first-mode 
disturbances  are  dominant  in  the  linear  and  nonlinear  stability  stages.  Two  breakdowns  at 
Mach  1.6  are  simulated  in  this  work:  fundamental  breakdown  and  oblique  breakdown. 

For  the  fundamental  breakdown,  we  use  the  following  parameters: 


Moo 

= 

1.6, 

Reh 

= 

15,000, 

Too 

300  A", 

Pr 

= 

0.72, 

0.75, 

= 

2.28, 

0.02, 

= 

0.001, 

C 

r~ 

0.25, 

^  TT  TT 

z 
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The  estimated  wave  angle  of  the  enforced  disturbance  \s  xj)  ^  60".  A  361  X  71  X  37  grid 
is  used,  with  the  wing  located  between  the  41st  and  321st  streamwise  grid  points.  The 
height  of  the  computational  domain  is  about  23/»,  and  the  chordlength/half-thickness  is 
60,  which  makes  the  Reynolds  number  (based  on  the  chordlength)  Rei  =  900, 000.  The 
blowing/suction  slot  is  located  between  the  61st  and  77th  grid  points.  For  the  time  stepping, 
one  forcing  period  (based  on  ui2d)  is  divided  into  600  time  steps.  It  takes  around  25  CRAY- 
C90  CPU  hours  and  25  MW  memory  for  our  code  to  run  14  forcing  periods. 

Figure  13  depicts  the  instantaneous  contour  plots  of  u',  v',  w',  p',  and  p'  after  14 
forcing  periods.  It  can  be  seen  that  even  though  we  imposed  only  a  pair  of  very  weak 
3-D  disturbances  (5%  the  magnitude  of  2-D’s),  those  disturbances  are  very  unstable  (see 
Figure  13  (c)  for  spanwise  component  of  perturbation  velocity).  The  process  of  transition 
is  still  the  same  as  that  in  subsonic  flow:  2D  and  3D  disturbance  development,  generation 
of  lambda  waves,  and  high  shear  layer  breakdown.  Not  as  we  expected  in  the  low  speed 
boundary  layer  flow,  the  spanwise  velocity  components  can  grow  even  in  the  favorable- 
pressure-gradient  regime  before  the  crest  point  of  the  wing,  where  the  2-D  disturbances 
decrease  in  this  regime.  A  quick  breakdown  occurs  after  the  disturbances  enter  the  adverse- 
pressure-gradient  regime,  and  the  disturbance  quantities  reach  the  saturated  level  very 
early.  This  property  can  be  observed  more  clearly  form  Figure  14,  in  which  we  give  the 
maximum  rms  value  of  the  perturbation  quantities  in  two  selected  (2,  y)  planes,  2  =  0  and 
2  =  ^.  The  quick  growth  of  the  disturbances  appears,  and  after  breakdown,  the  rms  of 
disturbances  keeps  at  a  relatively  stable  level  (saturated). 

The  instantaneous  contour  plots  of  the  streamwise  and  spanwise  perturbation  vortices 
as  well  as  the  perturbation  vorticity  magnitude  are  also  given  in  Figure  15,  which  provides 
a  better  understanding  of  the  fundamental  breakdown  structure  for  the  2D  airfoil  in  low 
supersonic  Mach  number. 

Another  case  simulated  in  this  work  is  the  “oblique  breakdown”.  For  this  case,  the 
same  parameters  as  for  the  simulation  of  fundamental  breakdown  are  used,  except  for 
the  following  differences:  At  the  blowing/suction  slot,  we  impose  only  one  pair  of  3D 
disturbance  waves  with  finite  amplitude  Sid  =  0.0075,  and  with  controlled  spanwise  wave 
number  p  =  1.59,  making  the  estimated  wave  angle  ip  =  45". 

As  in  the  fundamental  breakdown  case,  the  3D  disturbances  can  grow  in  the  favorable 
pressure  gradient  regime  before  they  reach  the  crest  point  of  the  wing.  This  is  different 
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from  incompressible  and  subsonic  flow  transition.  A  quick  breakdown  occurs  after  the 
disturbances  enter  the  adverse  pressure-gradient  regime.  Figure  16  gives  the  maximum 
rms  value  of  the  perturbation  quantities  in  two  selected  (x,y)  planes,  z  =  0  and  z  = 
Quick  growth  of  the  disturbances  is  observed.  By  comparing  Figure  14  and  Figure  16, 
we  And  that  the  disturbances  in  the  oblique  breakdown  case  increase  faster  than  those  in 
the  fundamental  case.  The  saturated  level  also  appears  earlier  than  in  the  fundamental 
breakdown  case.  Therefore,  we  can  consider  that  the  oblique  breakdown  is  a  major  factor 
for  boundary  layer  transition  at  low  supersonic  Mach  number. 

Figure  17  depicts  the  instantaneous  contour  plots  of  disturbance  quantities,  u',  u',  w\  p', 
and  p'.  Obviously,  the  structure  of  3D  disturbance  development  is  different  from  that  of  the 
fundamental  case.  No  lambda  wave  is  observed  in  this  case.  Instead,  a  “zig-zag”  type  of  the 
vortices  are  observed,  followed  by  a  branching  structure.  The  instantaneous  contour  plots 
of  the  streamwise  and  spanwise  perturbation  vortices  as  well  as  the  perturbation  vorticity 
magnitude  are  given  in  Figure  18,  which  provides  a  better  understanding  of  the  oblique 
breakdown  structure  for  the  2D  airfoil  in  low  supersonic  Mach  number.  It  can  also  be  found 
that  two  high-shear  layers  per  wavelength  in  the  spanwise  direction  appear,  organized  in 
the  streamwise  direction  in  a  staggered  pattern. 

9  Concluding  Remarks 

In  summary,  the  current  approach  shows  the  ability  of  DNS  to  simulate  the  whole  process  of 
flow  transition  around  2-D  airfoils  with  medium  Reynolds  number  («  10®)  for  both  subsonic 
and  low  supersonic  Mach  numbers.  It  demonstrates  that  realistic  numerical  simulations 
of  compressible  transition  phenomena  based  on  the  spatial  approach  are  already  feasible 
with  current  generation  of  supercomputers.  Although  the  computational  results  agree  well 
with  LST  and  qualitatively  agree  well  with  theory  and  experiment,  the  code  still  needs 
to  be  further  validated.  Once  validated,  this  simulation  can  be  extended  to  explore  new 
physical  phenomena  and  may  partiaJly  substitute  for  wind  tunnel  experiments,  especially 
for  high-speed  flows. 
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Figure  4.  (a)  Comparison  between  DNS  and  LST  of  the  disturbance  amplitude  of  the 
u  and  V  .  (b)  Comparison  of  DNS  and  LST  results  for  the  amplitude  eigenfunctions  of 
subsonic  flat  plate  boundary  layer.  Moo  =  0.5,  Re*  =  875,  a;  =  0.1. 


Figure  5.  (a)  Instantaneous  contour  plots  of  the  streamwise  perturbation  vorticity  on  the 
y~  =  0.18595,  0.48165,  0.85052  (x,2:)  planes  (from  bottom  to  top)  at  t  =  lOT  for  the  flat 
plate  fundamental  breakdown  case. 


Figure  5.  (b)  Instantaneous  contour  plots  of  the  spanwise  perturbation  vorticity  on  the 
y*  =  0.18595,  0.48165,  0.85052  (i,  z)  planes  (from  bottom  to  top)  at  t  =  lOT  for  the  flat  plate 
fundamental  breakdown  case.  Flow  parameter:  Moo  =  0.5,  Too  =  288iF,  Re*  -  875,  Pr  = 
0.72,  w  =  0.1,  e2i  =  0.02,  ezi  =  0.008,  ^  =  0.2.  Grids:  193  x  41  x  32. 
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Figure  6.  (a)  Instantaneous  contour  plots  of  the  streamwise  perturbation  vorticity  uix  on  the 
y'  =  0.1471,0.3768,0.6563  (x,2)  planes  (from  bottom  to  top)  at  t  =  IIT. 
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Figure  6.  (b)  Instantaneous  contour  plots  of  the  spanwise  perturbation  vorticity  Uz  on  the 
y*  =  0.1471,  0.3768,  0.6563  {x,z)  planes  (from  bottom  to  top)  at  t  =  llT  for  the  flat  plate 
subharmonic  breakdown  case.  Flow  parameter:  Moo  =  0.8,  Too  =  2SSK,  Re  =  1675,  Pr  = 
0.72,  U2d  =  0.15075,  £2<f  =  0.012,  Csd  =  0.004,  /?  =  0.22.  Grids:  265  x  45  x  42. 


Figure  7.  Grid  structure  for  the  subsonic  Joukowsky  airfoil  transition  simulation. 
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Streamwise  distribution  of  the  maximum  perturbation  rms  amplitudes  for  the 


Mach  0.5  Joukowsky  airfoil.  Reh  =  4000,  u)  =  0.8. 
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''■  Figure  9. 

--"^stantaneous  contour  plots  of 
y  /^a)  u',  (b)  u',  (c)  u;',  (d)  p',  (e)  p' 

y^--"  for  the  fundamental  breakdown  case 

of  a  Mach  0.5  airfoil  boundary  layer  on  the 
j  =  10  grid  surface  («  0.0641h)  at  t  =  lOT. 
Flow  parameter:  Moo  =  0.5,  Too  =  288K,  Reh  =  4000, 
Pr  =  0.72,  u)  =  0.8,  e^d  =  0.015,  e^d  =  0.01,  0  =  1.3. 
X  41  X  31,  chordlength/half-thickness=40  (Rei  =  160,000). 


Grids:  361  x  45  x  41,  chordlength/half-tliickness=50  {Rei  =  670,000). 
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Figure  11.  Compaxison  between  DNS  and  LST  of  the  disturbance  amplitude  of  the  second¬ 
mode  disturbance  u'  cind  v'  for  Mach  4.5  flat  plate  boimdary  layer.  Re*  =  10, 000,  uj  = 
2.0467. 
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Figure  12.  Amplitude  eigenfunctions  for  supersonic  flat  plate  transition:  comparison  of 
DNS  and  LST  results  for  the  second-mode  disturbance  at  Moo  =  4-5  and  Re*  =  10, 000. 
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Figure  13. 


Instantaneous  contour  plots 
(a)  u',  (b)  u',  (c)  w\  (d)  p',  (e)  p' 
^e-"^  for  the  fundamental  breakdown  case  of  a 
Mach  1.6  Joukowsky  airfoil  boundary  layer 
^  on  the  j  =  10  grid  surface  («  0.05/i)  at  t  =  12T. 

Flow  parameter:  M^o  =  1-6,  Too  =  SOO/f,  Ren  =  15000, 
Pr  =  0.72,  a;  =  0.75,  ead  =  0.02,  £3^  =  0.001,  (3  =  2.28. 
Grids:  361  x  71  x  37,  chordlength/half-thickness=60  (Rei  =  900,000). 


^o^io[(Maximum  rms  perturbation)  X  1000] 


Figure  14.  Streamwise  distribution  of  the  maximum  rms  disturbance  for  the  fundamental 
breakdown  case.  Moo  =  1-6,  F  =  50  x  10“®,  and  Rei  =  900,000.  Upper:  z  =  0.,  lower: 
2  =  .0689. 
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Figure  15.  Instantaneous  contour  plots  of  (a)  (b)  cj',  (c)  |a;'|  for  the  fundamental 

breakdown  case  of  a  Mach  1.6  Joukowsky  airfoil  boundary  layer  on  the  j  —  20  grid  surface 
(«  0.l2h  from  wall)  at  t  =  12T.  Flow  parameter:  Moo  =  1-6,  Too  —  300 iF,  Ren  =  15000, 
Pr  =  0.72,  Lo  =  0.75,  ejd  =  0.02,  £3d  =  0.001,  0  =  2.28.  Grids:  361  x  71  x  37, 
chordlength/half-thickness=60  (Rei  =  900,000). 
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Instantaneous  contour  plots 
of  (a)  u',  (b)  v',  (c)  w',  (d)  p',  (e)  p' 
for  the  oblique  breakdown  caise  of  a 
Mach  1.6  Joukowsky  airfoil  boundary  layer 
on  the  j  =  10  grid  surface  («  0.05/i)  at  i  =  12T. 
Flow  parameter:  Moo  =  1-6,  Too  =  300/f,  Reh  =  15000, 
Pr  =  0.72,  w  =  0.75,  e2d  =  0.0,  egd  =  0.0075,  0  =  1.59. 
Grids:  361  x  71  x  37,  chordlength/half-thickness=60  {Rei  =  900,000). 


Figure  18.  Instantaneous  contour  plots  of  (a)  (b)  a;',  (c)  |u;'|  for  the  oblique  breakdown 

case  of  a  Mach  1.6  Joukowsky  airfoil  boundary  layer  on  the  j  —  20  grid  surface  («  O.lh)  at 
t  =  12T.  Flow  parameter:  Moo  =  1-6,  Too  =  SOO/f,  Reh  =  15000,  Pr  =  0.72,  uj  =  0.75, 
e2d  =  0.0,  e^d  =  0.0075,  /?  =  1.59.  Grids:  361  x  71  x  37,  chordlength/half-thickness=60 
{Rei  =  900,000). 


